#########################################################################   
#                                 INFO                                  #   
#########################################################################  

# PROJECT: Exposure to State Violence and Substance Use
# AUTHOR: Alexandra Blackman; Elizabeth Nugent; Sarah Kammourh
# PURPOSE: analysis of drug use 2009-2014
# CREATED: October 2022
# INPUTS:  Panel SYPE 0914V3.dta


#########################################################################   
#                                 SETUP                                 #   
######################################################################### 

#.rs.restartR()
rm(list = ls()) 
setwd("~/Dropbox/Egypt_Alex Sarah Liz/")

######## PACKAGES ######## 

need <- c("doBy","foreign", "lubridate","extrafont", "car", "plyr","dplyr", "broom","xtable",
          "reshape", "reshape2", "openxlsx", "tidyr", "dplyr", "ggplot2", "ggthemes","tidyverse",
          "RColorBrewer", "rms", "lmtest", "sandwich", "msm","stargazer","texreg","memisc",
          "readr","purrr","tibble","stringr","forcats","jsonlite","gridExtra", "haven",
          "xml2","httr","rvest","DBI","hms","blob","magrittr","glue", "stringr",
          "readxl","multiwayvcov","readstata13","plm", "MASS", "dotwhisker") # packages needed
have <- need %in% rownames(installed.packages()) # installed packages
if(any(!have)) install.packages(need[!have]) # install missing packages
invisible(lapply(need, library, character.only=T)) # load needed packages

options(scipen=999)

#########################################################################   
#                 LOAD CLUSTERING, F-TEST FUNCTIONS                     #   
#########################################################################   

cluster.se <- function(model, cluster){
  coeftest(model, vcov=function(x) cluster.vcov(x, cluster))[,2]
}

f.test <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  ftest <- waldtest(model, vcov = vcovCL, test = "F")
  return(ftest)
} 

#########################################################################   
#                       LOAD DATA                                       #   
#########################################################################   

sype_panel <- read.dta13("Data/Panel SYPE 0914V3.dta")

###################################
#                                 #
#   DEPENDENT VARIABLES           #
#                                 #
###################################

######### TOBACCO ###########
# Q4207_14: Which statement do you think best describes your smoking behavior? (2014)
table(sype_panel$Q4207_14)
sype_panel$smokes_14 <- ifelse(sype_panel$Q4207_14=="Currently smoke occasionally" |
                                       sype_panel$Q4207_14=="Currently smoking" |
                                       sype_panel$Q4207_14=="Smoke tobacco products other than cigarettes daily",1,0)
table(sype_panel$smokes_14)
1327/(9576+1327) #12.2 percent have

# Q4207_09: Which statement do you think best describes your smoking behavior? (2009)
table(sype_panel$Q4207_09)
sype_panel$smokes_09 <- ifelse(sype_panel$Q4207_09=="i currently smoke occasionally" |
                                 sype_panel$Q4207_09=="currently smoking" |
                                 sype_panel$Q4207_09=="i smoke other tobacco products than cigarettes daily",1,0)
table(sype_panel$smokes_09)

# Q4203_4_14: Do your friends smoke? (2014)
table(sype_panel$Q4203_4_14)
sype_panel$friendssmoke_14 <- ifelse(sype_panel$Q4203_4_14=="All or most", 2, ifelse(sype_panel$Q4203_4_14 == "Few", 1, ifelse(sype_panel$Q4203_4_14 == "None", 0, NA)))
table(sype_panel$friendssmoke_14)

# Q4203_4_09: Do your friends smoke? (2009)
table(sype_panel$Q4203_4_09)
sype_panel$friendssmoke_09 <- ifelse(sype_panel$Q4203_4_09=="all/most", 2, ifelse(sype_panel$Q4203_4_09 == "few", 1, ifelse(sype_panel$Q4203_4_09 == "none", 0, NA)))
table(sype_panel$friendssmoke_09)


# Q4203_1_14: Does your father smoke? / Q4203_2_14: Does your mother smoke? (2014)
table(sype_panel$Q4203_1_14)
table(sype_panel$Q4203_2_14)
table(sype_panel$Q4203_3_14)

sype_panel$parentsmoke_14 <- ifelse(sype_panel$Q4203_1_14=="Yes" | sype_panel$Q4203_2_14=="Yes", 1, 0)
sype_panel$siblingsmoke_14 <- ifelse(sype_panel$Q4203_3_14=="All or most" | sype_panel$Q4203_3_14=="Few", 1, 0)

table(sype_panel$parentsmoke_14)
table(sype_panel$siblingsmoke_14)
sype_panel$familysmoke_14 <- ifelse(sype_panel$parentsmoke_14==1 | sype_panel$siblingsmoke_14==1, 1, 0)
table(sype_panel$familysmoke_14)

# Q4203_1_09: Does your father smoke? / Q4203_2_09: Does your mother smoke? (2009)
table(sype_panel$Q4203_1_09)
table(sype_panel$Q4203_2_09)
table(sype_panel$Q4203_3_09)

sype_panel$parentsmoke_09 <- ifelse(sype_panel$Q4203_1_09=="yes" | sype_panel$Q4203_2_09=="yes", 1, 0)
sype_panel$siblingsmoke_09 <- ifelse(sype_panel$Q4203_3_09=="all/most" | sype_panel$Q4203_3_09=="few", 1, 0)

table(sype_panel$parentsmoke_09)
table(sype_panel$siblingsmoke_09)
sype_panel$familysmoke_09 <- ifelse(sype_panel$parentsmoke_09==1 | sype_panel$siblingsmoke_09==1, 1, 0)
table(sype_panel$familysmoke_09)


######### ALCOHOL ###########
# Q4215_14: Have you tried any alcoholic drinks? (2014)
table(sype_panel$Q4215_14)
sype_panel$drinks_14 <- ifelse(sype_panel$Q4215_14=="No",0,1)
table(sype_panel$drinks_14)
241/(10662+241) #2.2 percent have or refuse to answer

# Q4215_09: Have you tried any alcoholic drinks? (2009)
table(sype_panel$Q4215_09)
sype_panel$drinks_09 <- ifelse(sype_panel$Q4215_09=="no",0,1)
table(sype_panel$drinks_09)

# Q4214_14: During the last year have any of your friends consumed any alcoholic drinks? (2014)
table(sype_panel$Q4214_14)
sype_panel$friendsdrink_14 <- ifelse(sype_panel$Q4214_14=="Yes",1,0)
table(sype_panel$friendsdrink_14)

# Q4214_09: During the last year have any of your friends consumed any alcoholic drinks? (2009)
table(sype_panel$Q4214_09)
sype_panel$friendsdrink_09 <- ifelse(sype_panel$Q4214_09=="yes",1,0)
table(sype_panel$friendsdrink_09)


######### DRUGS ###########
# Q4221_14: Have you experimented with any drugs before? (2014)
table(sype_panel$Q4221_14)
sype_panel$drugs_14 <- ifelse(sype_panel$Q4221_14=="No",0,1)
table(sype_panel$drugs_14)
194/(194+10709) #1.8 percent have

# Type of drug: Q4222_A_14 Q4222_B_14 Q4222_C_14 Q4222_D_14 Q4222_E_14 Q4222_F_14 Q4222_G_14 Q4222_H_14 Q4222_I_14 Q4222_J_14 Q4222_X_14 Q4222_K_14
table(sype_panel$Q4222_A_14) #Pills
table(sype_panel$Q4222_B_14) #Stimulants
table(sype_panel$Q4222_C_14) #Sedatives
table(sype_panel$Q4222_D_14) #Hashish
table(sype_panel$Q4222_E_14) #Heroine
table(sype_panel$Q4222_F_14) # Marijuana
table(sype_panel$Q4222_G_14) # Paint sniffing
table(sype_panel$Q4222_H_14) # Petrol Sniffing
table(sype_panel$Q4222_I_14) # Glue sniffing
table(sype_panel$Q4222_J_14) # Cocaine

# The vast majority of drug use is cannabis products -- hashish or banjo -- followed by pills (e.g., tramadol)

sype_panel$cannabis_14 <- ifelse(sype_panel$Q4222_D_14=="Yes" | sype_panel$Q4222_F_14=="Yes",1,0)
table(sype_panel$cannabis_14)

# Q4223_14: How often do you use drugs?
table(sype_panel$Q4223_14)

# Q4221_09: Have you experimented with any drugs before? (2009)
table(sype_panel$Q4221_09)
sype_panel$drugs_09 <- ifelse(sype_panel$Q4221_09=="no",0,1)
table(sype_panel$drugs_09)
111/(111+7650) #1.4 percent have

# Q4219_14: Have any of your friends experimented with drugs before?
table(sype_panel$Q4219_14)
sype_panel$friendsdrugs_14 <- ifelse(sype_panel$Q4219_14=="No",0,1)
table(sype_panel$friendsdrugs_14)

# Q4219_09: Have any of your friends experimented with drugs before?
table(sype_panel$Q4219_09)
sype_panel$friendsdrugs_09 <- ifelse(sype_panel$Q4219_09=="no",0,1)
table(sype_panel$friendsdrugs_09)

# Q4220_14: Have any of your family members experimented with drugs before?
table(sype_panel$Q4220_14)
sype_panel$familydrugs_14 <- ifelse(sype_panel$Q4220_14=="No",0,1)
table(sype_panel$familydrugs_14)

# Q4220_09: Have any of your family members experimented with drugs before?
table(sype_panel$Q4220_09)
sype_panel$familydrugs_09 <- ifelse(sype_panel$Q4220_09=="no",0,1)
table(sype_panel$familydrugs_09)


###################################
#                                 #
#   INDEPENDENT VARIABLE          #
#                                 #
###################################

## MAKE WITNESS VIOLENCE DURING REVOLUTION VARIABLE ##
# Think about the period since the start of the January 25th revolution. 
# Did you experience any of the following? Witness violence on TV
table(sype_panel$Q5417_A_14)
sype_panel$tvviolence_14 <- ifelse(sype_panel$Q5417_A_14=="Many times" |
                                      sype_panel$Q5417_A_14=="A few times" |
                                      sype_panel$Q5417_A_14=="Once",1,0)

table(sype_panel$tvviolence_14)

# Think about the period since the start of the January 25th revolution. 
# Did you experience any of the following? Witness violence on social media (Facebook, internet, etc)
table(sype_panel$Q5417_B_14)
sype_panel$smviolence_14 <- ifelse(sype_panel$Q5417_B_14=="Many times" |
                                   sype_panel$Q5417_B_14=="A few times" |
                                   sype_panel$Q5417_B_14=="Once",1,0)

table(sype_panel$smviolence_14)

# Think about the period since the start of the January 25th revolution. 
# Did you experience any of the following? Saw someone being injured in person
table(sype_panel$Q5417_C_14)

sype_panel$witness_injury_14 <- ifelse(sype_panel$Q5417_C_14=="Many times" |
                                      sype_panel$Q5417_C_14=="A few times" |
                                      sype_panel$Q5417_C_14=="Once",1,0)
table(sype_panel$witness_injury_14)

# Think about the period since the start of the January 25th revolution. 
# Did you experience any of the following? Saw someone killed in person
table(sype_panel$Q5417_D_14)

sype_panel$witness_death_14 <- ifelse(sype_panel$Q5417_D_14=="Many times" |
                              sype_panel$Q5417_D_14=="A few times" |
                              sype_panel$Q5417_D_14=="Once",1,0)
table(sype_panel$witness_death_14)

sype_panel$witnessviolence_14 <- ifelse(sype_panel$witness_injury_14==1 | sype_panel$witness_death_14==1,1,0)
table(sype_panel$witnessviolence_14)

# correlation matrix - exposure to different types of violence (2014)
tab_corr14 <- dplyr::select(sype_panel, witnessviolence_14, tvviolence_14, smviolence_14) %>% 
  na.omit() %>% 
  cor() %>% 
  xtable(caption= "Correlation Matrix: Exposure to Types of Violence (2014)")
digits(tab_corr14) <- 3
tab_corr14

sype_panel$violencetype_14 <- NA
sype_panel$violencetype_14 <- ifelse(sype_panel$tvviolence_14==1,"On TV",sype_panel$violencetype_14)
sype_panel$violencetype_14 <- ifelse(sype_panel$smviolence_14==1,"Social Media",sype_panel$violencetype_14)
sype_panel$violencetype_14 <- ifelse(sype_panel$witnessviolence_14==1,"In-Person",sype_panel$violencetype_14)
table(sype_panel$violencetype_14)

###################################
#                                 #
#   DEMOGRAPHIC CONTROLS          #
#                                 #
###################################

# Gender
table(sype_panel$sex_14)
sype_panel$resp_fem <- ifelse(sype_panel$sex_14=="female",1,0)

# Age in 2009
table(sype_panel$Q5_14)    #Age in 2009
table(sype_panel$Q6_14)    #Age in 2014
table(sype_panel$Q5_Y_14)  #Birth year

sype_panel$age_09 <- 2009 - sype_panel$Q5_Y_14
table(sype_panel$age_09)

# Age in 2014
sype_panel$age_14 <- 2014 - sype_panel$Q5_Y_14
table(sype_panel$age_14)

# Age^2 in 2009
sype_panel$age2_09 <- (sype_panel$age_09)^2
table(sype_panel$age2_09)

# Marital Status
table(sype_panel$Q7_14)
table(sype_panel$Q4111_14)
sype_panel$married_widowed_14 <- ifelse(sype_panel$Q7_14=="Contractually married" |
                                       sype_panel$Q7_14=="Married" |
                                       sype_panel$Q7_14=="Widowed",1,0)
table(sype_panel$married_widowed_14)

# Highest level of education attained
table(sype_panel$Q16_14)
table(sype_panel$Q17_14)

# Highest level of education attended (2014)
table(sype_panel$Q1101_14) #ever attend?
table(sype_panel$Q1105_14) 
sype_panel$school_14 <- ifelse(sype_panel$Q1101_14 == "Never been", "Never attended", ifelse(sype_panel$Q1105_14=="General primary" |
                                                                              sype_panel$Q1105_14=="Azhar primary" |
                                                                                sype_panel$Q1105_14=="General preparatory" |
                                                                                sype_panel$Q1105_14=="Azhar preparatory","Below secondary","Secondary+"))

sype_panel$school_14 <- factor(sype_panel$school_14, levels= c("Never attended","Below secondary","Secondary+"))
table(sype_panel$school_14)

# Highest level of education attended (2009)
table(sype_panel$q26_09)    #ever attend?
table(sype_panel$q27_09)    #highest level attended

sype_panel$school_09 <- ifelse(sype_panel$q26_09 == "never attended school", "Never attended", ifelse(sype_panel$q27_09=="primary" |
                                                                                               sype_panel$q27_09=="primary azhar" |
                                                                                               sype_panel$q27_09=="preparatory" |
                                                                                               sype_panel$q27_09=="preparatory azhar","Below secondary","Secondary+"))

sype_panel$school_09 <- factor(sype_panel$school_09, levels= c("Never attended","Below secondary","Secondary+"))
table(sype_panel$school_09)


# Employment Status
## Worked in last 7 days ##
table(sype_panel$Q2119_14)
table(sype_panel$Q2119_09)
sype_panel$employed_14 <- ifelse(sype_panel$Q2119_14=="yes",1,0)
sype_panel$employed_09 <- ifelse(sype_panel$Q2119_09=="yes",1,0)
table(sype_panel$employed_14)
table(sype_panel$employed_09)

## Unemployed ##
table(sype_panel$Q2122_14)
table(sype_panel$Q2122_09)
sype_panel$unemployed_14 <- ifelse(sype_panel$Q2119_14=="no" & sype_panel$Q2122_14=="yes",1,0)
sype_panel$unemployed_09 <- ifelse(sype_panel$Q2119_09=="no"& sype_panel$Q2122_09=="yes",1,0)
table(sype_panel$unemployed_14)
table(sype_panel$unemployed_09)

# wealth index quintiles 2014
table(sype_panel$wiq_14)
sype_panel$seswealth_14 <- ifelse(sype_panel$wiq_14=="Lowest",1,
                                   ifelse(sype_panel$wiq_14=="Second",2,
                                          ifelse(sype_panel$wiq_14=="Middle",3,
                                                 ifelse(sype_panel$wiq_14=="Fourth",4,5))))
table(sype_panel$seswealth_14)

# wealth index quintiles 2009
table(sype_panel$wiq_09)
sype_panel$seswealth_09 <- ifelse(sype_panel$wiq_09=="lowest",1,
                                   ifelse(sype_panel$wiq_09=="second",2,
                                          ifelse(sype_panel$wiq_09=="middle",3,
                                                 ifelse(sype_panel$wiq_09=="fourth",4,5))))
table(sype_panel$seswealth_09)

##Have you ever used the Internet? (2014 only) ##
table(sype_panel$Q5901_14)
sype_panel$internet_14 <- ifelse(sype_panel$Q5901_14=="Yes",1,0)

##Do you have internet at home?
table(sype_panel$Q5904_A_14)
table(sype_panel$Q5904_A_09)
sype_panel$internethome_14 <- ifelse(is.na(sype_panel$Q5904_A_14)==T, 0, ifelse(sype_panel$Q5904_A_14=="yes",1,0))
sype_panel$internethome_09 <- ifelse(is.na(sype_panel$Q5904_A_09)==T, 0, ifelse(sype_panel$Q5904_A_09=="yes",1,0))
table(sype_panel$internethome_14)
table(sype_panel$internethome_09)

## belong to a. community or associations/NGOs, b. sports/youth clubs, c. unions/syndicates, d. political party/movement/religion based but politically motivated e. unregistered, informal

sype_panel$groups_14 <- ifelse(sype_panel$Q5201_A_14=="yes"|sype_panel$Q5201_B_14=="yes"|sype_panel$Q5201_C_14=="yes"|sype_panel$Q5201_D_14=="yes"|sype_panel$Q5201_E_14=="yes",1,0)
table(sype_panel$groups_14)

## how satisfied are you with: the extent you feel like you belong in your community
table(sype_panel$Q5524_C_14)
sype_panel$belong <- ifelse(sype_panel$Q5524_C_14 == "Very satisfied", 5, 
          ifelse(sype_panel$Q5524_C_14 == "Satisfied", 4, 
          ifelse(sype_panel$Q5524_C_14 == "Neither satisfied nor dissatisfied", 3, 
                 ifelse(sype_panel$Q5524_C_14 == "Dissatisfied", 2, 
                        ifelse(sype_panel$Q5524_C_14 == "Very dissatisfied", 1, NA)))))
table(sype_panel$belong)

## how satisfied are you with: the social support you receive from your community
table(sype_panel$Q5524_E_14)

sype_panel$support <- ifelse(sype_panel$Q5524_E_14 == "Very satisfied", 5, 
                            ifelse(sype_panel$Q5524_E_14 == "Satisfied", 4, 
                                   ifelse(sype_panel$Q5524_E_14 == "Neither satisfied nor dissatisfied", 3, 
                                          ifelse(sype_panel$Q5524_E_14 == "Dissatisfied", 2, 
                                                 ifelse(sype_panel$Q5524_E_14 == "Very dissatisfied", 1, NA)))))
table(sype_panel$support)

# how many friends do you have? (2014)
table(sype_panel$Q5204_A_14)
sype_panel$friendsf_2014 <- ifelse(sype_panel$Q5204_A_14 < 98, sype_panel$Q5204_A_14, NA)

table(sype_panel$Q5204_B_14)
sype_panel$friendsm_2014 <- ifelse(sype_panel$Q5204_B_14 < 98, sype_panel$Q5204_B_14, NA)

sype_panel$nfriends_14 <- sype_panel$friendsf_2014 + sype_panel$friendsm_2014 
table(sype_panel$nfriends_14)

# how many friends do you have? (2009)
table(sype_panel$Q5204_A_09)
sype_panel$friendsf_2009 <- ifelse(sype_panel$Q5204_A_09 < 98, sype_panel$Q5204_A_09, NA)

table(sype_panel$Q5204_B_09)
sype_panel$friendsm_2009 <- ifelse(sype_panel$Q5204_B_09 < 98, sype_panel$Q5204_B_09, NA)

sype_panel$nfriends_09 <- sype_panel$friendsf_2009 + sype_panel$friendsm_2009 
table(sype_panel$nfriends_09)

###################################
#                                 #
#         DEPRESSION              #
#                                 #
###################################
#The DSM-5 outlines the following criterion to make a diagnosis of depression. 
#The individual must be experiencing five or more symptoms during the same 2-week period 
#and at least one of the symptoms should be either (1) depressed mood or (2) loss of interest or pleasure.
# Source: https://www.psychiatry.org/patients-families/depression/what-is-depression
# and https://www.psycom.net/depression-definition-dsm-5-diagnostic-criteria/

#Depressed mood most of the day, nearly every day.
## Depression: Do you feel unhappy?
table(sype_panel$Q4235_14)
table(sype_panel$Q4235_09)
sype_panel$depression1_14 <- ifelse(sype_panel$Q4235_14=="Yes",1,0)
sype_panel$depression1_09 <- ifelse(sype_panel$Q4235_09==1,1,0)
table(sype_panel$depression1_14)
table(sype_panel$depression1_09)

#Markedly diminished interest or pleasure in all, or almost all, activities most of the day, nearly every day.
## Depression: Do you find it difficult to enjoy your daily activities?
table(sype_panel$Q4237_14)
table(sype_panel$Q4237_09)
sype_panel$depression2_14 <- ifelse(sype_panel$Q4237_14=="Yes",1,0)
sype_panel$depression2_09 <- ifelse(sype_panel$Q4237_09==1,1,0)
table(sype_panel$depression2_14)
table(sype_panel$depression2_09)

#Significant weight loss when not dieting or weight gain, or decrease or increase in appetite nearly every day.
## Depression: Is your appetite poor?
table(sype_panel$Q4228_14)
table(sype_panel$Q4228_09)
sype_panel$depression3_14 <- ifelse(sype_panel$Q4228_14=="Yes",1,0)
sype_panel$depression3_09 <- ifelse(sype_panel$Q4228_09==1,1,0)
table(sype_panel$depression3_14)
table(sype_panel$depression3_09)

#A slowing down of thought and a reduction of physical movement (observable by others, not merely subjective feelings of restlessness or being slowed down).
## Depression: Do you have trouble thinking clearly?
table(sype_panel$Q4234_14)
table(sype_panel$Q4234_09)
sype_panel$depression4_14 <- ifelse(sype_panel$Q4234_14=="Yes",1,0)
sype_panel$depression4_09 <- ifelse(sype_panel$Q4234_09==1,1,0)
table(sype_panel$depression4_14)
table(sype_panel$depression4_09)

#Fatigue or loss of energy nearly every day.
## Depression: Do you feel tired all the time?
table(sype_panel$Q4244_14)
table(sype_panel$Q4244_09)
sype_panel$depression5_14 <- ifelse(sype_panel$Q4244_14=="Yes",1,0)
sype_panel$depression5_09 <- ifelse(sype_panel$Q4244_09==1,1,0)
table(sype_panel$depression5_14)
table(sype_panel$depression5_09)

#Feelings of worthlessness or excessive or inappropriate guilt nearly every day.
## Depression: Do you feel that you are a worthless person?
table(sype_panel$Q4242_14)
table(sype_panel$Q4242_09)
sype_panel$depression6_14 <- ifelse(sype_panel$Q4242_14=="Yes",1,0)
sype_panel$depression6_09 <- ifelse(sype_panel$Q4242_09==1,1,0)
table(sype_panel$depression6_14)
table(sype_panel$depression6_09)

#Diminished ability to think or concentrate, or indecisiveness, nearly every day.
## Depression: Do you find it difficult to make decisions?
table(sype_panel$Q4238_14)
table(sype_panel$Q4238_09)
sype_panel$depression7_14 <- ifelse(sype_panel$Q4238_14=="Yes",1,0)
sype_panel$depression7_09 <- ifelse(sype_panel$Q4238_09==1,1,0)
table(sype_panel$depression7_14)
table(sype_panel$depression7_09)

#Insomnia or hypersomnia
## Depression: Do you sleep badly?
table(sype_panel$Q4229_14)
table(sype_panel$Q4229_09)
sype_panel$depression8_14 <- ifelse(sype_panel$Q4229_14=="Yes",1,0)
sype_panel$depression8_09 <- ifelse(sype_panel$Q4229_09==1,1,0)
table(sype_panel$depression8_14)
table(sype_panel$depression8_09)

#Recurrent thoughts of death, recurrent suicidal ideation without a specific plan, or a suicide attempt or a specific plan for committing suicide.
## Depression: Has the thought of committing suicide been on your mind?
table(sype_panel$Q4243_14)
table(sype_panel$Q4243_09)
sype_panel$depression9_14 <- ifelse(sype_panel$Q4243_14=="Yes",1,0)
sype_panel$depression9_09 <- ifelse(sype_panel$Q4243_09==1,1,0)
table(sype_panel$depression9_14)
table(sype_panel$depression9_09)

## Depression Index

sype_panel$depressionindex_14 <- rowSums(sype_panel[, c("depression1_14","depression2_14",
                                                        "depression3_14","depression4_14",
                                                        "depression5_14","depression6_14",
                                                        "depression7_14","depression8_14",
                                                        "depression9_14")], na.rm = TRUE)

sype_panel$depressionindex_09 <- rowSums(sype_panel[, c("depression1_09","depression2_09",
                                                        "depression3_09","depression4_09",
                                                        "depression5_09","depression6_09",
                                                        "depression7_09","depression8_09",
                                                        "depression9_09")], na.rm = TRUE)

table(sype_panel$depressionindex_14)
table(sype_panel$depressionindex_09)
sype_panel$depression_14 <- ifelse(sype_panel$depressionindex_14>=5,1,0)
sype_panel$depression_09 <- ifelse(sype_panel$depressionindex_09>=5,1,0)
table(sype_panel$depression_14)
table(sype_panel$depression_09)

###################################
#                                 #
#         RELIGIOSITY             #
#                                 #
###################################

## Religion ##
table(sype_panel$q_105_09)
sype_panel$muslim <- NA
sype_panel$muslim[sype_panel$q_105_09 == "muslim"] <- 1
sype_panel$muslim[sype_panel$q_105_09 == "others"] <- 0
table(sype_panel$muslim)

## Q5801_A_14: Are you someone who prays daily? (1 = Never, 5 = Always) ##
table(sype_panel$Q5801_A_14)
sype_panel$pray_14 <- NA
sype_panel$pray_14[sype_panel$Q5801_A_14 == "Never"] <- 1
sype_panel$pray_14[sype_panel$Q5801_A_14 == "Rarely"] <- 2
sype_panel$pray_14[sype_panel$Q5801_A_14 == "Sometimes"] <- 3
sype_panel$pray_14[sype_panel$Q5801_A_14 == "Most of the time"] <- 4
sype_panel$pray_14[sype_panel$Q5801_A_14 == "Always"] <- 5

## Q5802: Regardless of whether you attend religious services or not,  ##
##        would you say that you are...1 = A very religious person, 3 = Not a religious person) ##
table(sype_panel$Q5802_14)
sype_panel$religious_14 <- NA
sype_panel$religious_14[sype_panel$Q5802_14 == "Not a religious person"] <- 1
sype_panel$religious_14[sype_panel$Q5802_14 == "A religious person"] <- 2
sype_panel$religious_14[sype_panel$Q5802_14 == "A very religious person"] <- 3

table(sype_panel$Q5802_09)
sype_panel$religious_09 <- NA
sype_panel$religious_09[sype_panel$Q5802_09 == "not a religious person"] <- 1
sype_panel$religious_09[sype_panel$Q5802_09 == "a religious person"] <- 2
sype_panel$religious_09[sype_panel$Q5802_09 == "a very religious person"] <- 3

## do you attend Friday prayers/Sunday prayers
table(sype_panel$Q5801_E_14)

sype_panel$communal_worship <- ifelse(sype_panel$Q5801_E_14 == "Never", 1, 
                             ifelse(sype_panel$Q5801_E_14 == "Rarely", 2, 
                                    ifelse(sype_panel$Q5801_E_14 == "Sometimes", 3, 
                                           ifelse(sype_panel$Q5801_E_14 == "Most of the time", 4, 
                                                  ifelse(sype_panel$Q5801_E_14 == "Always", 5, NA)))))
table(sype_panel$communal_worship)



###################################
#                                 #
#       OTHER BELIEFS             #
#                                 #
###################################

# Political interest: How often do you discuss politics with friends? (2009)
table(sype_panel$q_468_09)
sype_panel$polinterest_09 <- ifelse(sype_panel$q_468_09=="often", 3,
                                     ifelse(sype_panel$q_468_09=="sometimes", 2,
                                            ifelse(sype_panel$q_468_09=="rarely", 1, 0)))

# Political interest: How often do you discuss politics with friends? (2014)
table(sype_panel$Q5206_A_4_14)
sype_panel$polinterest_14 <- ifelse(sype_panel$Q5206_A_4_14=="Always" | sype_panel$Q5206_A_4_14=="Often", 3,
                                     ifelse(sype_panel$Q5206_A_4_14=="Sometimes", 2,
                                            ifelse(sype_panel$Q5206_A_4_14=="Rarely", 1, 0)))

table(sype_panel$polinterest_09)
table(sype_panel$polinterest_14)

# Political corruption: what is the level of corruption in public institutions? (2009/2014)
table(sype_panel$Q5517_09)
table(sype_panel$Q5517_14)

sype_panel$corruption_09 <- ifelse(sype_panel$Q5517_09==98, NA,sype_panel$Q5517_09)
sype_panel$corruption_14 <- ifelse(sype_panel$Q5517_14==98, NA,sype_panel$Q5517_14)

# Generalized Trust: would you say most people can be trusted? (2009/2014)
table(sype_panel$Q5302_09)
table(sype_panel$Q5302_14)

sype_panel$trust_09 <- ifelse(sype_panel$Q5302_09=="most people can be trusted",1,0)
sype_panel$trust_14 <- ifelse(sype_panel$Q5302_14=="Most people can be trusted",1,0)


# Political issues (2009)
#People need larger role in important governement decisions
sype_panel$people_role_09 <- ifelse(sype_panel$q_461_d_09=="most needed",1,0)

#Protecting freedom of speech (2009)
sype_panel$freedom_speech_09 <- ifelse(sype_panel$q_461_e_09=="most needed",1,0)

#Protecting political rights (2009)
sype_panel$political_rights_09 <- ifelse(sype_panel$q_461_j_09=="most needed",1,0)

sype_panel$demo_issues_09 <- sype_panel$people_role_09 + sype_panel$freedom_speech_09 + sype_panel$political_rights_09
table(sype_panel$demo_issues_09)

#Have you ever witnessed police violence? (2009)
# witnessing police violence in 2009 is used as baseline for 2014 violence variables (witness/tv/social media)
table(sype_panel$q_575_09)
sype_panel$policeviolence_09 <- ifelse(sype_panel$q_575_09=="yes",1,0)

sype_panel$witnessviolence_09 <- ifelse(sype_panel$q_575_09=="yes",1,0)

table(sype_panel$witnessviolence_09)

table(sype_panel$witnessviolence_14)

### set alternative 09 violence variables equal to police violence
sype_panel$tvviolence_09 <- sype_panel$witnessviolence_09
sype_panel$smviolence_09 <- sype_panel$witnessviolence_09

### Governorate Fixed Effects
table(sype_panel$GOV_14)
sype_panel$governorate <- factor(sype_panel$GOV_14)

### sampling unit for clustered standard errors
table(sype_panel$PSU_14)
sype_panel$psu <- factor(sype_panel$PSU_14)

sype_panel$ind_ID <- sype_panel$indid_14


###############################################################################
#                   EXPORT CLEAN SYPE DATA                                    #
#                                                                             #
###############################################################################
myvars <- c("witnessviolence_14","witnessviolence_09",
            "tvviolence_14","tvviolence_09",
            "smviolence_14","smviolence_09",
            "drugs_14","drugs_09",
            "smokes_14","smokes_09",
            "drinks_14","drinks_09",
            "friendssmoke_14","friendssmoke_09", 
            "familysmoke_14","familysmoke_09",
            "friendsdrink_14","friendsdrink_09",
            "friendsdrugs_14","friendsdrugs_09",
            "familydrugs_14","familydrugs_09",
            "age_14","age_09","school_14","school_09","employed_14","employed_09",
            "unemployed_14","unemployed_09","seswealth_14","seswealth_09","depression_14","depression_09",
            "religious_14","religious_09", "trust_14","trust_09", "nfriends_14","nfriends_09",
            "polinterest_14","polinterest_09","corruption_14","corruption_09", "violencetype_14",
            "ind_ID","resp_fem","muslim","married_widowed_14","policeviolence_09",
            "governorate","psu")

sype_vars <- sype_panel[myvars]
# write.xlsx(sype_vars, "Data/sype_vars.xlsx")

#########################################################################   
#               SUMMARY STATS ON MAIN IND AND DEP VARIABLES             #   
######################################################################### 

sumvars <- c("witnessviolence_14","witnessviolence_09",
              "drugs_14","drugs_09",
              "drinks_14","drinks_09",
              "smokes_14","smokes_09")

sype_sums <- sype_panel[sumvars]

# Make quick summary stats table
sum_stats <- sype_sums %>% 
  dplyr::summarize_all(mean, na.rm=T) %>% 
  xtable(caption= "Summary: Exposure to State Violence and Drug Use (2009 and 2014)")
digits(sum_stats) <- 3
sum_stats


#########################################################################   
#               DIFFERENCE IN MEANS BY EXPOSURE AND YEAR                #   
######################################################################### 

meanvars <- c("witnessviolence_14","witnessviolence_09",
               "drugs_14","drugs_09")

sype_means <- sype_panel[meanvars]

sype_means <- mutate(sype_means, witness_14 = ifelse(is.na(witnessviolence_14) == T,NA,
                                                     ifelse(witnessviolence_14 == 1, "Exposure", "No Exposure")))

sype_means <- mutate(sype_means, witness_09 = ifelse(is.na(witnessviolence_09) == T,NA,
                                                     ifelse(witnessviolence_09 == 1, "Exposure", "No Exposure")))

##### CREATE DIFFERENCE IN MEANS TABLE (2014)
table_means14 <- sype_means %>% 
  do(tidy(t.test(drugs_14~witnessviolence_14, data = .))) %>% 
  dplyr::rename(
    difference = estimate,
    mean_no_exp = estimate1,
    mean_exp = estimate2
  ) %>% 
  dplyr::select(-c(statistic, parameter, method, alternative,conf.low,conf.high))

table_means14 <- table_means14[c("mean_no_exp", "mean_exp", "difference","p.value")]

tab_means14 <-xtable(table_means14, 
                       caption= "Difference in Means: Exposure to State Violence (2014)")
digits(tab_means14) <- 3
tab_means14  


##### CREATE DIFFERENCE IN MEANS TABLE (2009)
table_means09 <- sype_means %>%
  do(tidy(t.test(drugs_09~witnessviolence_09, data = .))) %>%
  dplyr::rename(
    difference = estimate,
    mean_no_exp = estimate1,
    mean_exp = estimate2
  ) %>%
  dplyr::select(-c(statistic, parameter, method, alternative,conf.low,conf.high))

table_means09 <- table_means09[c("mean_no_exp", "mean_exp", "difference","p.value")]

tab_means09 <-xtable(table_means09,
                     caption= "Difference in Means: Exposure to State Violence (2009)")
digits(tab_means09) <- 3
tab_means09




###############################################################################
#                       MAKE LONG DATASET                                     #
#                                                                             #
###############################################################################

panelvars <- c("witnessviolence_14","witnessviolence_09",
               "tvviolence_14", "tvviolence_09", 
               "smviolence_14", "smviolence_09",
               "drugs_14","drugs_09",
               "smokes_14","smokes_09",
               "drinks_14","drinks_09",
               "friendssmoke_14","friendssmoke_09", 
               "familysmoke_14","familysmoke_09",
               "friendsdrink_14","friendsdrink_09",
               "friendsdrugs_14","friendsdrugs_09",
               "familydrugs_14","familydrugs_09",
               "age_14","age_09","school_14","school_09","employed_14","employed_09",
               "unemployed_14","unemployed_09","seswealth_14","seswealth_09","depression_14","depression_09",
               "religious_14","religious_09", "trust_14","trust_09", "nfriends_14","nfriends_09",
               "polinterest_14","polinterest_09","corruption_14","corruption_09", 
               "ind_ID","resp_fem","muslim","governorate","psu")

sype_wide <- sype_panel[panelvars]

#Remove witness violence & drug NAs
panel_nocontrols <- sype_wide[complete.cases(sype_wide[ , 1:4]),]
sype_long_nocontrols <-reshape(panel_nocontrols, varying=c(1:44), direction="long", idvar="ind_ID", sep="_", timevar="year")

#Remove witness violence & drug NAs & all controls
panel_controls <- sype_wide[complete.cases(sype_wide[ , 1:4]),]
panel_controls <- panel_controls[complete.cases(panel_controls[ , 15:40]),]
sype_long_controls <-reshape(panel_controls, varying=c(1:44), direction="long", idvar="ind_ID", sep="_", timevar="year")

# head(sype_long_controls)
# head(select(sype_panel, ind_ID2, witnessviolence_14,witnessviolence_09,tvviolence_14,tvviolence_09,smviolence_14,smviolence_09))
# head(select(p_controls, ind_ID2, witnessviolence,tvviolence,smviolence))

###############################################################################
#                        PANEL ANALYSIS FOR DRUGS                             #
#                                                                             #
###############################################################################
# No controls
head(sype_long_nocontrols)
sype_long_nocontrols$years <- ifelse(sype_long_nocontrols$year==14,"2014","2009")
sype_long_nocontrols$years <- factor(sype_long_nocontrols$years, levels= c("2009","2014"))

class(sype_long_nocontrols$ind_ID)
sype_long_nocontrols$ind_ID2 <- as.factor(sype_long_nocontrols$ind_ID)

col_order <- c("ind_ID2", "years","ind_ID", "resp_fem", "muslim", "governorate", "psu", "year", 
               "witnessviolence", "tvviolence", "smviolence", "drugs", "smokes", "drinks", "friendssmoke", "familysmoke", 
               "friendsdrink", "friendsdrugs","familydrugs", "age","school", "employed", "unemployed",
               "seswealth", "depression", "religious", "trust", "nfriends","polinterest", "corruption")
sype_long_nocontrols <- sype_long_nocontrols[, col_order]
sype_long_nocontrols

p_nocontrols <- pdata.frame(sype_long_nocontrols, index=c("ind_ID2","years"), drop.index=F, row.names=TRUE)


# Controls
head(sype_long_controls)
sype_long_controls$years <- ifelse(sype_long_controls$year==14,"2014","2009")
sype_long_controls$years <- factor(sype_long_controls$years, levels= c("2009","2014"))

class(sype_long_controls$ind_ID)
sype_long_controls$ind_ID2 <- as.factor(sype_long_controls$ind_ID)

col_order <- c("ind_ID2", "years","ind_ID", "resp_fem", "muslim", "governorate", "psu", "year", 
               "witnessviolence", "tvviolence", "smviolence", "drugs", "smokes", "drinks", "friendssmoke", "familysmoke", 
               "friendsdrink", "friendsdrugs","familydrugs", "age","school", "employed", "unemployed",
               "seswealth", "depression", "religious", "trust", "nfriends","polinterest", "corruption")
sype_long_controls <- sype_long_controls[, col_order]
sype_long_controls

p_controls <- pdata.frame(sype_long_controls, index=c("ind_ID2","years"), drop.index=F, row.names=TRUE)

### PANEL ANALYSIS FOR DRUGS

# Create list of model names
models.panel <- list()

# Run models
# two-way fixed effects panel models (individual and year fixed effects)
models.panel[[1]] <- plm(drugs ~ witnessviolence, data=p_nocontrols, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel[[2]] <- plm(drugs ~ witnessviolence + friendsdrugs + familydrugs + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")

names(models.panel) <- c("Drug Use", "Drug Use")

# Cluster-robust Huber/White standard errors by ind_ID2
models.panel.cl <- list()
models.panel.cl[[1]] <- coeftest(models.panel[[1]], vcov=vcovHC(models.panel[[1]], cluster="group"))[,2]
models.panel.cl[[2]] <- coeftest(models.panel[[2]], vcov=vcovHC(models.panel[[2]], cluster="group"))[,2]

stargazer(models.panel, 
          out = "Tables/plm_drugs.tex",
          title = "Panel Analysis: Exposure to Violence and Drug Use",
          label = "plm_drugs",
          type = "text", 
          se = models.panel.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("ind_ID"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.panel),
          dep.var.labels = c("Probability of Reported Use"),
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by respondent.",
                    "Model 2 also controls for time-varying measures of: ",
                    "friends' and family members' usage, employment, education,",
                    "SES, depression, religiosity, generalized trust, and number",
                    "of friends."))


###############################################################################
#                 PANEL ANALYSIS FOR ALCOHOL AND TOBACCO                      #
#                                                                             #
###############################################################################

#Remove witness violence & drug NAs
panel_nocontrols_app <- sype_wide[complete.cases(sype_wide[ , 1:8]),]
sype_long_nocontrols_app <-reshape(panel_nocontrols_app, varying=c(1:44), direction="long", idvar="ind_ID", sep="_", timevar="year")


#Remove witness violence & drug NAs & all controls
panel_controls_app <- sype_wide[complete.cases(sype_wide[ , 1:40]),]
sype_long_controls_app <-reshape(panel_controls_app, varying=c(1:44), direction="long", idvar="ind_ID", sep="_", timevar="year")

# No controls
head(sype_long_nocontrols_app)
sype_long_nocontrols_app$years <- ifelse(sype_long_nocontrols_app$year==14,"2014","2009")
sype_long_nocontrols_app$years <- factor(sype_long_nocontrols_app$years, levels= c("2009","2014"))

class(sype_long_nocontrols_app$ind_ID)
sype_long_nocontrols_app$ind_ID2 <- as.factor(sype_long_nocontrols_app$ind_ID)

col_order <- c("ind_ID2", "years","ind_ID", "resp_fem", "muslim", "governorate", "psu", "year", 
               "witnessviolence", "tvviolence", "smviolence", "drugs", "smokes", "drinks", "friendssmoke", "familysmoke", 
               "friendsdrink", "friendsdrugs","familydrugs", "age","school", "employed", "unemployed",
               "seswealth", "depression", "religious", "trust", "nfriends","polinterest", "corruption")
sype_long_nocontrols_app <- sype_long_nocontrols_app[, col_order]
sype_long_nocontrols_app

p_nocontrols_app <- pdata.frame(sype_long_nocontrols_app, index=c("ind_ID2","years"), drop.index=F, row.names=TRUE)


# Controls
head(sype_long_controls_app)
sype_long_controls_app$years <- ifelse(sype_long_controls_app$year==14,"2014","2009")
sype_long_controls_app$years <- factor(sype_long_controls_app$years, levels= c("2009","2014"))

class(sype_long_controls_app$ind_ID)
sype_long_controls_app$ind_ID2 <- as.factor(sype_long_controls_app$ind_ID)

col_order <- c("ind_ID2", "years","ind_ID", "resp_fem", "muslim", "governorate", "psu", "year", 
               "witnessviolence", "tvviolence", "smviolence", "drugs", "smokes", "drinks", "friendssmoke", "familysmoke", 
               "friendsdrink", "friendsdrugs","familydrugs", "age","school", "employed", "unemployed",
               "seswealth", "depression", "religious", "trust", "nfriends","polinterest", "corruption")
sype_long_controls_app <- sype_long_controls_app[, col_order]
sype_long_controls_app

p_controls_app <- pdata.frame(sype_long_controls_app, index=c("ind_ID2","years"), drop.index=F, row.names=TRUE)


### PANEL ANALSIS FOR ALCOHOL AND TOBACCO 

# Create list of model names
models.panel.app <- list()

# Run models
models.panel.app[[1]] <- plm(drinks ~ witnessviolence, data=p_nocontrols_app, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel.app[[2]] <- plm(drinks ~ witnessviolence + friendsdrink + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls_app, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel.app[[3]] <- plm(smokes ~ witnessviolence, data=p_nocontrols_app, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel.app[[4]] <- plm(smokes ~ witnessviolence + friendssmoke + familysmoke + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls_app, index=c("ind_ID2","years"), effect="twoways", model="within")

names(models.panel.app) <- c("Alcohol Use", "Alcohol Use", "Tobacco Use", "Tobacco Use")

# Cluster-robust Huber/White standard errors by ind_ID2
models.panel.app.cl <- list()
models.panel.app.cl[[1]] <- coeftest(models.panel.app[[1]], vcov=vcovHC(models.panel.app[[1]], cluster="group"))[,2]
models.panel.app.cl[[2]] <- coeftest(models.panel.app[[2]], vcov=vcovHC(models.panel.app[[2]], cluster="group"))[,2]
models.panel.app.cl[[3]] <- coeftest(models.panel.app[[3]], vcov=vcovHC(models.panel.app[[3]], cluster="group"))[,2]
models.panel.app.cl[[4]] <- coeftest(models.panel.app[[4]], vcov=vcovHC(models.panel.app[[4]], cluster="group"))[,2]

stargazer(models.panel.app, 
          out = "Tables/plm_alcohol_tobacco.tex",
          title = "Panel Analysis: Exposure to Violence and Substance Use (Alcohol and Tobacco)",
          label = "plm_alcohol_tobacco",
          type = "text", 
          se = models.panel.app.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("ind_ID"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.panel.app),
          dep.var.labels = c("Probability of Reported Use","Probability of Reported Use"),
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by respondent. Models 2 and 4 also",
                    "control for time-varying measures of: friends' and family members' usage,",
                    "employment, education, SES, depression, religiosity, generalized trust, and",
                    "number of friends."))

## COMBINE TABLES FOR PANEL ANALYSIS DRUG / ALCOHOL / TOBACCO USE

models.panel.main <- append(models.panel, models.panel.app)
models.panel.main.cl <- append(models.panel.cl, models.panel.app.cl)

cov.labels <- c("Witnessed Violence", "Friend Drug Use", "Family Drug Use",
                "Friend Alcohol Use", "Friend Tobacco Use", "Family Tobacco Use",
                "Employed", "Unemployed", "Attended School (Completed below secondary)",
                "Attended School (Completed secondary+)", "Wealth Quintile", "Depression",
                "Religiosity", "Generalized Trust", "Number of Friends")

stargazer(models.panel.main, 
          out = "Tables/plm_main.tex",
          title = "Panel Analysis: Exposure to Violence and Drug Use",
          label = "plm_main",
          type = "text", 
          se = models.panel.main.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("ind_ID"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.panel.main),
          dep.var.labels = rep("Probability of Reported Use", 3),
          covariate.labels = cov.labels,
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by respondent.",
                    "Models 2, 4, and 6 also control for time-varying measures of: ",
                    "friends' and family members' usage, employment, education,",
                    "SES, depression, religiosity, generalized trust, and number",
                    "of friends."))


###############################################################################
#                COEFFICIENT PLOT OF THE MARGINAL EFFECT                      #
#                                                                             #
###############################################################################

broom::tidy(models.panel.main[[2]])[1,]

fit1_df <- broom::tidy(models.panel.main[[2]])[1,] %>% 
  mutate(term = "Drugs",
         model = "m1",
         std.error = models.panel.main.cl[[2]][[1]])
fit2_df <- broom::tidy(models.panel.main[[4]])[1,] %>% 
  mutate(term = "Alcohol",
         model = "m1",
         std.error = models.panel.main.cl[[4]][[1]])
fit3_df <- broom::tidy(models.panel.main[[6]])[1,] %>% 
  mutate(term = "Tobacco",
         model = "m1",
         std.error = models.panel.main.cl[[6]][[1]])

dwplot(rbind(fit1_df,fit2_df,fit3_df))+
  scale_x_continuous(labels = scales::percent_format())+
  labs(x="Estimated Increase in Use from Exposure to Violence",
       caption="Point estimates with 95% confidence intervals.")+
  geom_vline(xintercept = 0,
             colour = "grey60",
             linetype = 2)+
  theme_minimal()+
  theme(legend.position = "none")
ggsave("Figures/coefplot.png", bg="white", width = 5, height = 4)




###############################################################################
#                PANEL ANALYSIS W/ DIFFERENT EXPOSURE                         #
#                         TO VIOLENCE                                         #
###############################################################################

# Create list of model names
models.panel <- list()

models.panel[[1]] <- plm(drugs ~ tvviolence + friendsdrugs + familydrugs + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")
models.panel[[2]] <- plm(drugs ~ smviolence + friendsdrugs + familydrugs + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel[[3]] <- plm(drinks ~ tvviolence + friendsdrink + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")
models.panel[[4]] <- plm(drinks ~ smviolence + friendsdrink + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")

models.panel[[5]] <- plm(smokes ~ tvviolence + friendssmoke + familysmoke + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")
models.panel[[6]] <- plm(smokes ~ smviolence + friendssmoke + familysmoke + employed + unemployed + school +
                           seswealth + depression + religious + trust + nfriends, 
                         data=p_controls, index=c("ind_ID2","years"), effect="twoways", model="within")

names(models.panel) <- c("Drug Use","Drug Use","Alcohol Use","Alcohol Use","Tobacco Use","Tobacco Use")

# Cluster-robust Huber/White standard errors by ind_ID2
models.panel.cl <- list()
models.panel.cl[[1]] <- coeftest(models.panel[[1]], vcov=vcovHC(models.panel[[1]], cluster="group"))[,2]
models.panel.cl[[2]] <- coeftest(models.panel[[2]], vcov=vcovHC(models.panel[[2]], cluster="group"))[,2]
models.panel.cl[[3]] <- coeftest(models.panel[[3]], vcov=vcovHC(models.panel[[3]], cluster="group"))[,2]
models.panel.cl[[4]] <- coeftest(models.panel[[4]], vcov=vcovHC(models.panel[[4]], cluster="group"))[,2]
models.panel.cl[[5]] <- coeftest(models.panel[[5]], vcov=vcovHC(models.panel[[5]], cluster="group"))[,2]
models.panel.cl[[6]] <- coeftest(models.panel[[6]], vcov=vcovHC(models.panel[[6]], cluster="group"))[,2]

cov.labels.type <- c("Saw Violence on TV", "Saw Violence on Social Media",
                     "Friend Drug Use", "Family Drug Use",
                "Friend Alcohol Use", "Friend Tobacco Use", "Family Tobacco Use",
                "Employed", "Unemployed", "Attended School (Completed below secondary)",
                "Attended School (Completed secondary+)", "Wealth Quintile", "Depression",
                "Religiosity", "Generalized Trust", "Number of Friends")


stargazer(models.panel, 
          out = "Tables/plm_type.tex",
          title = "Panel Analysis: Exposure to Violence (Different Types of Exposure)",
          label = "plm_type",
          type = "text", 
          se = models.panel.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("ind_ID"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.panel),
          dep.var.labels = rep("Probability of Reported Use", 3),
          covariate.labels = cov.labels.type,
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by respondent.",
                    "All models also control for time-varying measures of: ",
                    "friends' and family members' usage, employment, education,",
                    "SES, depression, religiosity, generalized trust, and number",
                    "of friends."))

###############################################################################
#                CROSS-SECTIONAL ANALYSIS W/ 2009 CONTROLS                    #
#                                                                             #
###############################################################################

### CROSS-SECTIONAL FOR DRUGS 

# Create list of model names
models.crosssec <- list()

# Run models
models.crosssec[[1]] <- lm(drugs_14~witnessviolence_14 + resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec[[2]] <- lm(drugs_14~witnessviolence_14 + drugs_09 + friendsdrugs_09 + familydrugs_09 + 
                             employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                             trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                             governorate, data=sype_vars)

names(models.crosssec) <- c("Drug Use", "Drug Use")

# Cluster standard errors by psu
models.crosssec.cl <- list()
models.crosssec.cl[[1]] <- cluster.se(models.crosssec[[1]], sype_vars$psu)
models.crosssec.cl[[2]] <- cluster.se(models.crosssec[[2]], sype_vars$psu)


stargazer(models.crosssec, 
          out = "Tables/crosssec_main.tex",
          title = "Cross-Sectional Analysis: Exposure to Violence and Drug Use",
          label = "crosssec_main",
          type = "text", 
          se = models.crosssec.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("governorate","resp_fem", "age_14","employed_09", "unemployed_09", "school_09", "seswealth_09", 
                   "depression_09", "polinterest_09", "trust_09", "nfriends_09","religious_09", "muslim"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.crosssec),
          dep.var.labels = c("Probability of Reported Use"),
          #covariate.labels = c(), 
          add.lines = list(c("Sex/Age/Religion Controls", rep("Y", length(models.crosssec))),
                           c("Additional 2009 Controls", "N","Y"),
                           c("Governorate FE", rep("Y", length(models.crosssec)))), 
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by sampling unit.",
                    "All models control for respondent sex, age, and religion.",
                    "Model 2 also controls for 2009 measures of: employment, ",
                    "education, SES, depression, generalized trust, number of",
                    "friends, and religiosity."))



### CROSS-SECTIONAL FOR ALCOHOL AND TOBACCO 

# Create list of model names
models.crosssec.app <- list()

# Run models
models.crosssec.app[[1]] <- lm(drinks_14~witnessviolence_14 + resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec.app[[2]] <- lm(drinks_14~witnessviolence_14 + drinks_09 + friendsdrink_09 + 
                             employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                             trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                             governorate, data=sype_vars)

models.crosssec.app[[3]] <- lm(smokes_14~witnessviolence_14 + resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec.app[[4]] <- lm(smokes_14~witnessviolence_14 + smokes_09 + friendssmoke_09 + familysmoke_09 + 
                             employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                             trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                             governorate, data=sype_vars)

names(models.crosssec.app) <- c("Alcohol Use", "Alcohol Use", "Tobacco Use", "Tobacco Use")

# Cluster standard errors by psu
models.crosssec.app.cl <- list()
models.crosssec.app.cl[[1]] <- cluster.se(models.crosssec.app[[1]], sype_vars$psu)
models.crosssec.app.cl[[2]] <- cluster.se(models.crosssec.app[[2]], sype_vars$psu)
models.crosssec.app.cl[[3]] <- cluster.se(models.crosssec.app[[3]], sype_vars$psu)
models.crosssec.app.cl[[4]] <- cluster.se(models.crosssec.app[[4]], sype_vars$psu)


stargazer(models.crosssec.app, 
          out = "Tables/crosssec_alcohol_tobacco.tex",
          title = "Cross-Sectional Analysis: Exposure to Violence and Substance Use (Alcohol and Tobacco)",
          label = "crosssec_alcohol_tobacco",
          type = "text", 
          se = models.crosssec.app.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("governorate","resp_fem", "age_14","employed_09", "unemployed_09", "school_09", "seswealth_09", 
                   "depression_09", "polinterest_09", "trust_09", "nfriends_09","religious_09", "muslim"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.crosssec.app),
          dep.var.labels = c("Probability of Reported Use","Probability of Reported Use"),
          #covariate.labels = c(), 
          add.lines = list(c("Sex/Age/Religion Controls", rep("Y", length(models.crosssec.app))),
                           c("Additional 2009 Controls", "N","Y","N","Y"),
                           c("Governorate FE", rep("Y", length(models.crosssec.app)))), 
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by sampling unit. All models control",
                    "for respondent sex, age, and religion. Models 2, 4 and 6 also control for 2009 measures",
                    "of: employment, education, SES, depression, generalized trust, number of friends,",
                    "and religiosity."))


## COMBINE TABLES FOR CROSS SECTIONAL ANALYSIS DRUG / ALCOHOL / TOBACCO USE

models.cross.main <- append(models.crosssec, models.crosssec.app)
models.cross.main.cl <- append(models.crosssec.cl, models.crosssec.app.cl)

cov.labels <- c("Witnessed Violence (2014)", "Drug Use (2009)", 
                "Friend Drug Use (2009)", "Family Drug Use (2009)",
                "Alcohol Use (2009)", "Friend Alcohol Use (2009)", 
                "Tobacco Use (2009)","Friend Tobacco Use (2009)", "Family Tobacco Use (2009)",
                "Witnessed Violence (2009)")

stargazer(models.cross.main, 
          out = "Tables/crosssec_all.tex",
          title = "Exposure to Violence and Substance Use",
          label = "crosssec_all",
          type = "text", 
          se = models.cross.main.cl, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("governorate","resp_fem", "age_14","employed_09", "unemployed_09", "school_09", "seswealth_09", 
                   "depression_09", "polinterest_09", "trust_09", "nfriends_09","religious_09", "muslim"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.cross.main),
          dep.var.labels = rep("Probability of Reported Use", 3),
          covariate.labels = cov.labels,
          add.lines = list(c("Sex/Age/Religion Controls", rep("Y", length(models.cross.main))),
                           c("Additional 2009 Controls", "N","Y","N","Y","N","Y"),
                           c("Governorate FE", rep("Y", length(models.cross.main)))), 
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by sampling unit. All models control",
                    "for respondent sex, age, and religion. Models 2, 4, and 6 also control for 2009 measures",
                    "of: employment, education, SES, depression, generalized trust, number of friends,",
                    "and religiosity."))


###############################################################################
#                CROSS-SECTIONAL ANALYSIS W/ 2014 CONTROLS                    #
#                                                                             #
###############################################################################

# Create list of model names
models.crosssec14 <- list()

# Run models
models.crosssec14[[1]] <- lm(drugs_14~witnessviolence_14 + friendsdrugs_14 + familydrugs_14 + 
                             employed_14 + unemployed_14 + school_14 + seswealth_14 + depression_14 + 
                             trust_14 + nfriends_14 + religious_14 + muslim + resp_fem + age_14 + 
                             governorate, data=sype_vars)

models.crosssec14[[2]] <- lm(drinks_14~witnessviolence_14 + friendsdrink_14 +
                               employed_14 + unemployed_14 + school_14 + seswealth_14 + depression_14 + 
                               trust_14 + nfriends_14 + religious_14 + muslim + resp_fem + age_14 + 
                               governorate, data=sype_vars)

models.crosssec14[[3]] <- lm(smokes_14~witnessviolence_14 + friendssmoke_14 + familysmoke_14 +
                               employed_14 + unemployed_14 + school_14 + seswealth_14 + depression_14 + 
                               trust_14 + nfriends_14 + religious_14 + muslim + resp_fem + age_14 + 
                               governorate, data=sype_vars)

names(models.crosssec14) <- c("Drug Use", "Alcohol Use", "Tobacco Use")

# Cluster standard errors by psu
models.crosssec.cl14 <- list()
models.crosssec.cl14[[1]] <- cluster.se(models.crosssec14[[1]], sype_vars$psu)
models.crosssec.cl14[[2]] <- cluster.se(models.crosssec14[[2]], sype_vars$psu)
models.crosssec.cl14[[3]] <- cluster.se(models.crosssec14[[3]], sype_vars$psu)


stargazer(models.crosssec14, 
          out = "Tables/crosssec_2014.tex",
          title = "Cross-Sectional Analysis: Exposure to Violence and Substance Use (2014 Controls)",
          label = "crosssec_2014",
          type = "text", 
          se = models.crosssec.cl14, 
          omit.stat = c("f","adj.rsq","ser"),  
          omit = c("governorate"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.crosssec14),
          dep.var.labels = c("","Probability of Reported Use",""),
          #covariate.labels = c(), 
          add.lines = list(c("Governorate FE", rep("Y", length(models.crosssec14)))), 
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by sampling unit."))


###############################################################################
#                CROSS-SECTIONAL ANALYSIS W/ DIFFERENT EXPOSURE               #
#                         TO VIOLENCE                                         #
###############################################################################
# Create list of model names
models.crosssec.type <- list()

# Run models
models.crosssec.type[[1]] <- lm(drugs_14~witnessviolence_14 + tvviolence_14 + smviolence_14 + 
                             resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec.type[[2]] <- lm(drugs_14~witnessviolence_14 + tvviolence_14 + smviolence_14 + 
                                  drugs_09 + friendsdrugs_09 + familydrugs_09 + 
                                  employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                                  trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                                  governorate, data=sype_vars)

models.crosssec.type[[3]] <- lm(drinks_14~witnessviolence_14 + tvviolence_14 + smviolence_14 + 
                             resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec.type[[4]] <- lm(drinks_14~witnessviolence_14 + tvviolence_14 + smviolence_14 +
                                  drinks_09 + friendsdrink_09 +
                                  employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                                  trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                                  governorate, data=sype_vars)

models.crosssec.type[[5]] <- lm(smokes_14~witnessviolence_14 + tvviolence_14 + smviolence_14 + 
                             resp_fem + age_14 + muslim + governorate, data=sype_vars)

models.crosssec.type[[6]] <- lm(smokes_14~witnessviolence_14 + tvviolence_14 + smviolence_14 +
                                  smokes_09 + friendssmoke_09 + familysmoke_09 + 
                                  employed_09 + unemployed_09 + school_09 + seswealth_09 + depression_09 + 
                                  trust_09 + nfriends_09 + policeviolence_09 + religious_09 + muslim + resp_fem + age_14 + 
                                  governorate, data=sype_vars)

names(models.crosssec.type) <- c("Drug Use", "Drug Use", "Alcohol Use", 
                                 "Alcohol Use", "Tobacco Use", "Tobacco Use")

# Cluster standard errors by psu
models.crosssec.cl.type <- list()
models.crosssec.cl.type[[1]] <- cluster.se(models.crosssec.type[[1]], sype_vars$psu)
models.crosssec.cl.type[[2]] <- cluster.se(models.crosssec.type[[2]], sype_vars$psu)
models.crosssec.cl.type[[3]] <- cluster.se(models.crosssec.type[[3]], sype_vars$psu)
models.crosssec.cl.type[[4]] <- cluster.se(models.crosssec.type[[4]], sype_vars$psu)
models.crosssec.cl.type[[5]] <- cluster.se(models.crosssec.type[[5]], sype_vars$psu)
models.crosssec.cl.type[[6]] <- cluster.se(models.crosssec.type[[6]], sype_vars$psu)


stargazer(models.crosssec.type, 
          out = "Tables/crosssec_type.tex",
          title = "Cross-Sectional Analysis: Exposure to Violence and Substance Use (Different Types of Exposure)",
          label = "crosssec_type",
          type = "text", 
          se = models.crosssec.cl.type, 
          omit.stat = c("f","adj.rsq","ser"), 
          omit = c("governorate","resp_fem", "age_14","employed_09", "unemployed_09", "school_09", "seswealth_09", 
                   "depression_09", "polinterest_09", "trust_09", "nfriends_09", "religious_09", "muslim"),
          font.size = "scriptsize",
          multicolumn = T,
          column.labels = names(models.crosssec.type),
          dep.var.labels = c("","Probability of Reported Use",""),
          #covariate.labels = c(), 
          add.lines = list(c("Sex/Age/Religion Controls", rep("Y", length(models.crosssec.type))),
                           c("Additional 2009 Controls", "N","Y","N","Y","N","Y"),
                           c("Governorate FE", rep("Y", length(models.crosssec.type)))), 
          notes.align = "l",
          notes = c("OLS models with standard errors clustered by sampling unit. All models control for respondent sex, age, and religion.",
                    "Models 2, 4, and 6 also control for 2009 measures of: employment, education, SES, depression, generalized trust,",
                    "number of friends, and religiosity."))



stop("Analysis complete.")




